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Abstract 

We work out the non-equilibrium steady state properties of a harmonic lattice 
which is connected to heat reservoirs at different temperatures. The heat reservoirs 

■ are themselves modeled as harmonic systems. Our approach is to write quantum 
Q . Langevin equations for the system and solve these to obtain steady state properties 
CJ I such as currents and other second moments involving the position and momentum 

operators. The resulting expressions will be seen to be similar in form to results 
. obtained for electronic transport using the non-equilibrium Green's function formal- 

ism. As an application of the formalism we discuss heat conduction in a harmonic 
\^ ' chain connected to self-consistent reservoirs. We obtain a temperature dependent 

■ thermal conductivity which, in the high-temperature classical limit, reproduces the 
, exact result on this model obtained recently by Bonetto, Lebowitz and Lukkarinen. 
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1 Introduction 



The harmonic crystal is one of the the simplest model that one learns in solid state physics 
and it is known to reproduce, for example, some of the correct experimental features for 
the specific heat of an insulating solid. The harmonic approximation basically involves 
expanding the full atomic potential of the solid about its minimum (which one assumes 
is a crystal) and keeping terms up to second order. If one transforms to normal mode 
coordinates then the harmonic crystal can be viewed as a collection of noninteracting 
phonons. While many equilibrium properties can be understood satisfactorily within the 
harmonic approximation, transport properties (heat conduction) of the harmonic lattice 
are anomalous because of absence of interactions between the phonons. 

In discussions of transport properties (see for example Ref. PJ) of an insulating solid 
it is usual to first think of the pure harmonic crystal. One can then introduce various 
mechanisms for scattering of the phonons, two important ones being impurity scattering 
and phonon-phonon scattering. Impurity scattering can arise because of randomness in 
the masses of the particles or in the spring constants. In this case the phonons still do 
not interact with each other. One can think of the phonons of the original pure crystal 
getting elastically scattered by impurities. Alternatively, since the system is still harmonic, 
one can think of heat transmission by the new normal modes of the disordered system. 
The second mechanism for scattering is through phonon-phonon interactions and this 
occurs if we include the higher order nonlinear terms {i.e beyond quadratic order) of the 
interparticle potential. Phonon-phonon interactions are usually classified into those which 
conserve momentum and those which do not (Umklapp processes). 

One important question that has attracted a lot of attention is: what are the necessary 
and sufficient conditions for the validity of Fourier's law of heat conduction [21 [3] ? We 
recall that Fourier's law states that for a solid with a spatially varying temperature field 
T(x) inside it, the local heat current density J at a point x is given by: 



where k defines the thermal conductivity of the solid and is expected to be an intrinsic 
property of the material. Fourier's law is a phenomenological law which is expected to 
be true in the hydrodynamic linear response regime. However till now there does not 
exist any purely mechanical model (without external potentials) in which a first principle 
demonstration of the validity of Fourier's law, either numerically or analytically, has been 
achieved. 

Infact it is now pretty much clear that in one dimensional momentum conserving 
systems Fourier's law is not valid and one cannot define a system-size independent thermal 
conductivity for these systems [H El E] ( Note that the momentum one is referring to 
here is the total real momentum, and not the crystal momentum. This is conserved 
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if there are no external potentials). Apart from a large number of numerical studies 
various theoretical approaches have been used for different classes of systems to arrive 
at this conclusion regarding non-validity of Fourier's law in one dimensions. In the case 
of interacting systems (nonlinear) the anomalous behaviour of thermal conductivity has 
been understood within the Green-Kubo formalism and has been related to long-time tails 
in the current-current auto-correlation functions E]- For non- interacting (harmonic) 
disordered systems heat transport occurs through independent phonon modes and the 
main contribution comes from low frequency extended modes [7]. One finds a situation 
similar to that of electronic transport in mesoscopic systems [H [9] where it is important 
to include the reservoirs in discussing transport. In this case it is more natural to talk 
of heat conductance of a system and one finds that this depends on the details of the 
reservoirs. 

For a three dimensional solid, based on kinetic theory and Boltzmann equation ap- 
proaches, the expectation is that the thermal conductivity should be finite (and Fourier's 
law valid) at high temperatures where Umklapp processes dominate. At low tempera- 
tures Umklapp processes are rare and non-Umklapp ones cannot lead to a finite thermal 
conductivity. However it is expected that, if in addition to these non-Umklapp processes, 
one also includes impurity scattering , then one should get a finite conductivity [T]. Un- 
fortunately there is no clear proof of any of these physically motivated expectations. One 
possible approach to this problem is through a rigorous study of the phonon Boltzmann 
equation (see for example [TO]). 

In this paper we discuss a formalism for transport in harmonic lattices based on the 
Langevin equation approach. This approach was first used to study heat conduction in 
a one-dimensional ordered harmonic lattice Subsequently this approach was used 

to study heat conduction in disordered harmonic lattices in one [IH |T3l |7j and two di- 
mensions IHj. The quantum mechanical case has also been studied by several authors 
[13 [m [IZl dSl [ISl EDI EI] using an open system description either through quantum 
Langevin equations or through density matrices. The open system description for har- 
monic systems closely resembles the Landauer formalism used for electron transport and 
Ref. [20] gives a derivation of the basic Landauer results using the Langevin equation 
approach, both for electrons and phonons. Another rigorous approach to studying elec- 
tron transport in mesoscopic systems is the non-equilibrium Green's function formalism 
(NEGF)[8j and Ref. [22] shows how this can be derived, for non-interacting electrons 
modeled by tight-binding Hamiltonians, using a quantum Langevin equation approach. 
We note that the Landauer formalism has also been discussed in the context of wave 
propagation in disordered media [231 l2l] • The problems of heat conduction in disordered 
harmonic lattices and wave propagation in disordered media are closely related. It is ex- 
pected that some of the work in the latter area, for example on localization, will be useful 
in the heat conduction context. 
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In the present paper we show how the quantum Langevin equations method for har- 
monic lattices leads to NEGF-like expressions for phonon transport. For simplicity we 
restrict ourselves to harmonic Hamiltonians with scalar displacement variables at each 
lattice site. It is straightforward to extend the calculations to the case of vector dis- 
placements. The basic steps in the calculation are: (i) one thinks of the full system as 
consisting of the sample we are interested in (henceforth called wire) as well as the reser- 
voirs at different temperatures which are connected to the wire, (ii) the wire and the 
reservoir Hamiltonians are taken to be harmonic, (iii) we eliminate the reservoir degrees 
of freedom and this leads to Langevin equations of motion for the wire variables, (iv) 
the hnear Langevin equations are solved and steady state properties such as expectation 
values of the current are found. Finally (v) the solution is written in a form where one 
can identify the usual phonon Green's functions commonly used in solid state physics. 
This leads to the identification with results from the NEGF formalism. We note that 
Landauer-like results for phonons have been proposed earlier [25l |26] and a recent paper 
j27j derives NEGF results for phonon transport using the Keldysh approach. 

The paper is organized as follows. In Sec. [2] we present the basic model and derive 
the quantum Langevin equations describing the wire. In Sec. Elwe solve the equations of 
motion to obtain the stationary long-time solution which is used in Sec. IHto obtain formal 
expressions for various steady state quantities such as the heat current. In Sec.[5]we discuss 
an application of the present approach to studying heat conduction in a harmonic wire 
with self-consistent reservoirs at all lattice points. Finally we conclude with a discussion 
in Sec. El 

2 Quantum Langevin equations 

We consider a harmonic system which consists of a wire (denoted by W) coupled to 
reservoirs which are also described by harmonic interactions. In most of our discussions 
we consider the case of two reservoirs, labeled as L (for left) and R (right) , which are at 
two different temperatures. It is easy to generalize the case where there are more than 
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two reservoirs. The Hamiltonian of the entire system of wire and reservoirs is taken to be 



n = ^X^MX + ^X^$X (2.1) 

= Hw + nL + HR + VL + VR 

1 • • 1 

where Hw = -X^MwXw + -X'^^w^w , 

Hl = ^XIMlXl + ^XI^lXl , 

Hr = ^XIMrXr + ^XI<I>rXr, 

Vl = X^VlXi, Vr = X^VrXr , 

where M, M]y, Ml, Mr are real diagonal matrices representing masses of the particles 
in the entire system, wire, left, and right reservoirs respectively. The quadratic potential 
energies are given by the real symmetric matrices $, ^l, while Vl and Vr 

denote the interaction between the wire and the two reservoirs. The column vectors 
X, X]y, Xl, Xr are Heisenberg operators which correspond to particle displacements, 
assumed to be scalars, about some equilibrium configuration. Thus X = {Xi, X2, ...X^^}'^ 
where Xr denotes the position operator of the r*^ particle and Ng denotes the number of 
points in the entire system. Also X = M~^ P where Pr denotes the momentum operator, 
with {Xr, Pr} satisfying the usual commutation relations [Xr,-Ps] = ihSrs- 
The Heisenberg equations of motion for the system are: 

MwXw = -^wXw - VlXl - VrXr , (2.2) 
and the equations of motion for the two reservoirs are 

MlXl = -^lXl - V[Xw , (2.3) 
MrXr = -^rXr-V^Xw (2.4) 

We solve these equations by considering them as linear inhomogeneous equations. Thus 
for the left reservoir the general solution to Eq. 02.31) is (for t > to): 

XL{t) = fl{t-to)MLXL{to)+gt{t~to)MLXL{to) 

- [ dt' gl{t~t')VlXw{t') , (2.5) 

With /+(t) = ULCosin^ml 9{t), gl{t) = uJ^^p^Uld{t) , 
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where 9{t) is the Heaviside function and Ul, are the normal mode eigenvector and 
eigenvalue matrices respectively and which satisfy the equations: 

UI<^lUl = ^1, UIMlUl = I. 
Similarly, for the right reservoir we obtain 

Xnit) = f+{t-to)MRXR{to)+gUt-to)MRXR{to) 

- f dt' gl{t-t')ViXw{t') . (2.6) 

We plug these solutions back into the equation of motion for the system to get 

MwXw = -^wXw + Vl + [ dt' Vl gt{t-t')V[ Xw{t') 

J to 

+ m + [ dt' Vr g+it - t') Xwit') , (2.7) 

J to 

where 

VL = -VL[ft{t-h)MLXL{h) + gt{t-h)MLXL{h)] 

VR = -VR[f^{t-to)MRXR{to) + gl{t-h)MRXR{to)]. (2.8) 

This equation has the form of a quantum Langevin equation. The properties of the noise 
terms r]L and r]R are determined using the condition that, at time the two isolated 
reservoirs are described by equilibrium phonon distribution functions. At time to the left 
reservoir is in equilibrium at temperature and the population of the normal modes (of 
the isolated left reservoir) is given by the distribution function /^(a;, Tl) — l/fe'^/'^^^^ — 1]. 
The equilibrium correlations are then given by: 

{XL{to)Xl{t,)) = U,^coM^)Ul . 

2Ul 2kBiL 

{X,{t,)Xl{to)) = Ul^coM^WI 
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Using these we can determine the correlations of the noise terms in Eq. fl2.8p . Thus we 
get for the left reservoir noise correlations: 



{VL{t)vl{t')) = VlU^ 

h 



cosVtjit — t')—— coth (— — ^) 



-i sin ^lit — t') 



and a similar expression for the right reservoir. 



UlVl , (2.9) 



3 Stationary solution of the equations of motion 

Now let us take the limits of infinite reservoir sizes and let to ~^ — oo. We can then solve 
Eq. (12.71) by taking Fourier transforms. Thus defining the Fourier transforms 

1 

Xw{to) = — dt Xv^(t)e*"* , 



27r 



— oo 
oo 

iujt 



1 f°° 

2vr J-oo 



oo 



gtA^) = / dtgtAt)e'-\ (3.1) 

J —CO 

we get from Eq. (12.71) 

{-uj^ Mw + ^w)Xw{uj) = [J:1{uj) + J:^{uj)] Xw{oo) + fiL{uj) + f]R{io) 

where S+(a;) = V^(?+(^)Vf, i:^^) = VRg^iuj)V^ . (3.2) 

The noise correlations can be obtained from Eq. (12.91) and we get (for the left reservoir): 

IT 

= 6iu + u') Tl{u) -[1 + Mu, Tl)] (3.3) 
vr 

where Tl{uj) = /m[S+(cj)] (3.4) 

which is a fluctuation-dissipation relation. This also leads to the more commonly used 
correlation: 

i( fiL{u;)f^l{uj') + U^')vl{^) ) = 6{uj + uj') T^u;) ^ coth(-^). (3.5) 

2 ZTT IKbI L 



Similar relations hold for the noise from the right reservoir. We then get the following 
stationary solution to the equations of motion: 



/oo 
dujXw{oi!)e~ 
-oo 



■iuit 



with Xwiuj) = G^{uj) [riLiuj) + f]niuj)] , (3.6) 
where Gtr = tt—. , , . (3.7) 

The identification of G'^{uj) as a phonon Green function, with S^^(u;) as effective self 
energy terms, is the main step that enables a comparison of results derived by the quantum 
Langevin approach with those obtained from the NEGF method. In App. \M we show 
explicitly how this identification is made. 

For the reservoirs we get, from Eqs. (I2.5tl2.6p . 

~VlXl{uj) = f,L{uj) + i:lXw{uj) , 

-VrXr{uj) = 77^(cu) + S+Xw.(cu) . (3.8) 

4 Steady state properties 

Current: The simplest way to evaluate the steady state current is to evaluate the follow- 
ing expectation value for left-to-right current: 

duo / duo' e-''(-+-')*^cu {XI{uo)VlXl{uo')) , 

-co J — oo 

which is just the rate at which the left reservoir does work on the wire. Using the solution 
in Eq. fl3.6ti3.8p we get 



J = - duo duo' e-(-+-')*za; ( ( fjl{uo) + f^{uj) ) G+ ^(^) 

J ~oo J —oo 

X [ u^') + s+(^') (^0 ( u^') + mW) )]). (4.1) 

Now consider that part of J, say J^, which depends only on T/^. Clearly this is: 

/oo roo 
duo / du' e-*('^+^')* I uo 
oo J —oo 

xTr [ G^\uj) Etiuo') GU^') ( m{^')fni^) ) ] 

r°° T huo 

= - I duoTrl G^^iuo) ^U-uo)G^{-uo)Tr{uo)] — [1 + Muo,Tn> 

J-oo 71" 



Using the identities Gy^/^ = G^^, cj) = G^^{uj) and taking the real part of above 

equation we obtain, after some simphfications, 

/OO f- 
duj Tr[ G+ (^) TlH G^iuj) TRiuj) ] — [1 + MuJ,Tn)] . 
-OO 

Including the contribution from the terms involving T^, and noting that the current has 
to vanish for Tl = Tr, it is clear that the net current will be given by 

duj Tr[ G+ (cu) Tl{u;)G^{u) Tn{uj)] — [f{u, n) - f{u, Tr)] . (4.2) 

-OO 

This expression for current can be seen to be of identical form as the NEGF expression 
for electron current (see for example [8t [28l [29] ). 

Two point correlation functions: We can also easily compute expectation values 
of various correlations. Thus the velocity-velocity correlations are given by 

K = {XwX^) 



TX 



+ G+,{uj)TR{uj)G^{uj)nw{l + h{uj, Tn)) 
duj- [G^{uj)TLiuj)G:^{uj)^coth{- ^ 



+ G+ (c.)r«(c.)G^(c.)^ coth(-^) ] , (4.3) 

where the last line is easily obtained after writing K = {K + K*)/2. We see that for 
Tl = Tpt this reduces to the equilibrium result of Eq. flB.3l) provided that there are no 
hound states. 

Similarly the position-position and position-velocity correlations are given by: 
P = {XwX'^) 

= r du^ [ GUu^)rL{u)G^{u) coth(-^) 

J-oo 2.11 2kBlL 

+ GU^)Tr{u)G^{u) cotH- ^ ^ ' 



'2kBTR' 

C = {X]yX^) 



/OO • f- 

du- [G^{UJ)VL{UJ)G^{UJ)—C0ih{ 
■OO 71" 2 



2kBTL 



GU^)Vr{u)G-^{u)^ ^oth(-^) ] . (4.4) 
z zkb-Lr 
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The correlation functions K and P can be used to define the local energy density which 
can in turn be used to define the temperature profile in the non-equilibrium steady state 
of the wire. Also we note that the correlations C give the local heat current density. In 
the next section we will find that it is sometimes more convenient to evaluate the total 
steady state current from this expression rather than the one in Eq. (14. 2p . 

Classical limits: The classical limit is obtained by taking the high temperature limit 
so that huj/ksT —>■ 0. Then we obtain the following expressions for the various steady 
state properties computed in the last section. The current is given by 

J = ^^"^ ~ ^""^ r du Tr[ G^u^) Tl{u;)G^{u;) Tn{u;)] , (4.5) 

TT J-oo 

while other correlation functions are given by: 



K = Mi: / dco CO G^{uj)Tl{uj)G^{u;) 

-oo 

fOO 

did uj G^{uj)V R{uj)Gy^{uj) 

'OO 



IT 



ksTL r , 1 



./-oo ^ 



vr 



duj - G^{uj)Tl{uj)G^{uj) 

CO 

I duj - G^{uj)TR{uj)Gyy{uj) 

J~oo ^ 



oo 
oo 



C = I du G^{u)Vl{uj)G^{u) 



TT 
TT 



-oo 

/oo 
du G^{Lu)TRiLu)G^{Lu) . (4.6) 
-oo 



For one dimensional wires these lead to [20] expressions for current and temperature used 
in earlier studies of heat conduction in disordered harmonic chains p!2| [T3l [7]. 

5 An application: one-dimensional harmonic crystal 
with self-consistent reservoirs 

As an application of the Langevin equation-Green functions formalism we consider the 
problem of heat transport in a harmonic chain with each site connected to self-consistent 
heat reservoirs. The classical version of this model was first studied by [501 [5T] who 
introduced the self-consistent reservoirs as a simple scattering mechanism which might 
ensure local equilibration and the validity of Fourier's law. This model was recently solved 
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exactly by Bonetto et al [33] who proved local equilibration and validity of Fourier's law 
and obtained an expression for the thermal conductivity of the wire. They also showed 
that the temperature profile in the wire was linear. The quantum version of the problem 
was also studied by Visscher and Rich [32] who analyzed the limiting case of weak coupling 
to the self-consistent reservoirs. We will show here how the present formalism can be used 
to obtain results in the quantum-mechanical case. The classical results of Bonetto et 
al are obtained as the high temperature limit while the quantum mechanical results of 
Vischer and Rich are obtained in the weak coupling limit. 

In this model one considers a one-dimensional harmonic wire described by the Hamil- 
tonian 



N N+1 2 

1=1 1=1 



2 

5 



= ^X^MwXw + ^X^^wXw (5.1) 

where the wire particles are denoted as X"^ = {xi,X2, ■■■Xisf} and we have chosen the 
boundary conditions Xq = xjy+i = 0. All the particles are connected to heat reservoirs 
which are taken to be Ohmic. The coupling strength to the reservoirs is controlled by 
the dissipation constant 7. The temperatures of the first and last reservoirs are fixed 
and taken to be Ti = and T/v = Tr. For other particles, i.e I = 2,3...{N — 1), the 
temperature of the attached reservoir T; is fixed self-consistently in such a way that the 
net current flowing into any of the reservoirs / = 2,3...(A^ — 1) vanishes. The Langevin 
equations of motion for the particles on the wire are: 

mxi = -mujl{2xi - xi-i - xi^i) - mujlxi - '^xi + r]i I = 1,2...N , (5.2) 

where the noise-noise correlation is easier to express in frequency domain and given by 

J( rii{uj)ri^{uj') + r]i{uj')r]miuj) ) = coth( ) S{uj + uj') Sim ■ (5.3) 

2 Zn ZKb-li 

From the equations of motion it is clear that the l^^ particle is connected to a bath with a 
self energy matrix S^(ct;) whose only non vanishing element is = i'yuj. Generalizing 

Eq. (14. 2 p to the case of multiple baths we find that the heat current from the reservoir 
into the wire is given by: 

Ji = J2 ^wi^) Vi{uj)G^{u) r^(u;)]— [f{u,Ti) - f{u,Tm)] , (5.4) 

m=l ^ 

where = [ -uj'^ Mw + ^w-^ S+(cu) , F; = /m[S+] . 
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Using the form of Ti we then get: 



To find the temperature profile we need to solve the N — 2 nonlinear equations J/ = 
for / = 2,3...A^ — 1 with Ti = and T/v = Tr. To proceed we consider the linear 
response regime with the applied temperature difference AT = — Tr << T where 
T = (Tl + Tji)/2. In that case we expand the phonon distribution functions f{u!,Ti) 
about the mean temperature T and get the following simpler expressions for the currents 

Ji = 7 / — of — 2^ I [Gw{^)]lm I [Tl - Tra) • (5.5) 

m=l 



We write = Z^^/{muj1) where Z is a tridiagonal matrix with offdiagonal elements 
equal to —1 and diagonal elements are all equal to z = 2 + uJq/ujI — 1^1 — i^oj/ {mujl). 
It is then easy to find its inverse using the formula Z^^ = Di^i^i Dm+i,N / Di^m^ where 
Dij is the determinant of the sub-matrix of Z beginning with the ith row and col- 
umn and ending with the jth row and column. The determinant is given by Dij = 
sinh [(j — i + 2) a]/ sinh(Q;), where e° = z/2 ± [(z/2)^ — 1]^/^ (any of the two roots can be 
taken). For points far from the boundaries of the wire (/ = yN where y = 0(1), 1 — y = 
0(1)) we then find that 

^—a\l—m\ 

" 2mcu2sinha ' ^^'^^ 

where we choose the root a such that = Re[a] > 0. In Ref. [33], where the classical 
version of the present model was studied, it was shown that in the limit N ^ oo the 
temperature profile obtained by solving the self-consistent equations has the linear form 

Ti = TL + jf^iTn~TL) . (5.7) 

From the form of in Eq. (15.61) we see at once that, for any point / in the bulk of the 
wire, the zero-current condition J/ = is satisfied since ^m=-oo(^ ~ m)|e~"l'~'"lp = 0. 
For points which are within distance 0(1) from the boundaries the temperature profile 
deviates from the linear form. 

To find the net left-right current in the wire we could use the formula for Ji = —Jn 
given in Eq. 05.41] . However we notice that use of this formula requires us to know the 
accurate form of for points m close to the boundaries since these terms contribute 
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significantly to the sum in Eq. fl5.4p . We instead use a different expression for tlie current. 
We evaluate the left-right current Ji^i+i on the bond connecting sites / and (/ + !). Using 
Eq. (14.41) and in the linear response limit, making expansions about T, we get 

T 2/ • \ [°° , f hw y ^ hw 

Ji^ i+i = mu^ {xixi+i) = / duj uj [ cosech 



IT \2kBT J 2kBT' 



oo 

X 

r?i=— oo 



The current value is independent of / and we choose to evaluate it at a value of I in the 
bulk of the wire. In that case terms in the sum above which contain with m close to 
the boundaries are exponentially small (~ e~°^) and so do not contribute. Hence we can 
use the linear temperature profile for and the form of in Eq. (15.61) to evaluate the 
heat current. We get: 

7 [°° , ^ f huj Y 2, ^ ^ 

J = —- ^ — : / duj -— — — — — cosech 



muj'^m J IsinhaP \2kBT J 2kBT 



^ kBTn 



where the additional terms in the summation over m are exponentially small contributions. 
Also as before we choose the root a{u) such that Re[a] > 0. Using the notation ur^u) = 
Re[a], ai{uj) = Im[a] we get, after some algebra, the following expression for the thermal 
conductivity k = JN/AT (obtained in the large N limit): 

K = — / du 5 — — cosech^ — — — . 5.8 

IGmu^in sinh an \2kBT J zft^T ysmha smha*/ 

In the high temperature limit we get (/ia;/2fcB7')^cosech^(/ia;/2A;BT) 1. This gives the 
classical result for the thermal conductivity. In this limit a change of variables from uj to 
aj leads to the following result for the thermal conductivity: 

2fcBmw2(2 + z/2) sin2(a/) 
Kci — / daj 



77r Jo (2 + z/2)2 _4cos2(aj) 

kBmuj'^ 



7 (2 + Z/2 + [z/2(4 + Z/2)]l/2) ' 

where z/ = ujq/uc- This agrees with the result obtained in Ref. [33] . 



(5.9) 
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An interesting limiting case is the case of weak coupling to the reservoirs (7 — 0). In 
this case Eq. f l5.8l) gives: 

— ^ — / da/ sin^a/ cosech^(— — ^) , (5.10) 
ksT J 477r 2kBT 

where cu^ = tUg + 2co'^[l — cos(a/)] . 

This agrees with the result obtained in [32]. The temperature profile obtained by them 
differs from the linear form obtained by us and also by [33] and is incorrect. They 
nevertheless obtain the correct thermal conductivity, presumably because their derivation 
only uses the temperature profile close to the centre of the chain where it is again linear. 
In the low temperature limit, Eq. (IS.lOp gives ~ Q-fi^^o/ksT ^ q ^^^^^ 
Ky^c ~ T for ujq = 0. As noted in [32] the expression for thermal conductivity (in the 
weak scattering limit) is consistent with a simple relaxation-time form for the thermal 
conductivity. The temperature dependence of then simply follows the temperature 
dependence of the specific heat of the 1-dimensional chain. 

Finally we examine the general case where the coupling constant has a finite value. 
In Figs. fni2l) we plot the thermal conductivity as a function of temperature for two sets 
of parameter: (i) •y/^muJc) = 0.2, u = 0.5 and (ii) •y/^mUc) = 0.2, = 0.0. The insets 
in the two figures show the low-temperature behaviour. As before, the low-temperature 
behaviour depends on whether or not there is an onsite potential. However we see that the 
form of the low-temperature behaviour is very different from the case of weak coupling. 
For small T it is easy to pull out the temperature dependence of the integral in Eq. (15. 8p 
and we find that k ~ for z/ 7^ and k, ~ for z/ = 0, in agreement with the 
numerical result shown in Figs. (11121) . 

6 Discussion 

In this paper, using the quantum Langevin equation approach, we have derived NEGF- 
like expressions for the heat current in a harmonic lattice connected to external reservoirs 
at different temperatures. We note that unlike other approaches such as the Green-Kubo 
formalism and Boltzmann equation approach, the Langevin equation approach explicitly 
includes the reservoirs. The Langevin equation is also physically appealing since it gives a 
nice picture of the reservoirs as sources of noise and dissipation. Also just as the Landauer 
formalism and NEGF have been extremely useful in understanding electron transport in 
mesoscopic systems it is likely that a similar description will be useful for the case of heat 
transport in insulating nanotubes, nanowires, etc. We note here that the single-channel 
Landauer results follow from NEGF if one considers one- dimensional reservoirs [2U] and 
have been useful in interpreting experimental results [M] . 
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We think that the similarity, between heat conduction studies in harmonic systems 
and electron transport in noninteracting wires, is an interesting and useful point to note. 
The two areas have developed quite independently using different theoretical tools. In 
the former case most of the earlier studies were done on classical systems using either a 
Langevin or a Fokker Planck description. More recent studies on quantum systems have 
used either a quantum Langevin or a density matrix approach. On the other hand, for 
the electron case, which is inherently quantum-mechanical, the most popular and useful 
approach has been the Landauer and the NEGF formalism. As we have demonstrated, in 
this paper for the phonon case, and in Ref [22] for the electron case, the NEGF results 
can be easily derived using the Langevin approach at least in the noninteracting case. 

As an application of the NEGF results we have studied the problem of heat conduction 
in a one dimensional harmonic chain connected to self-consistent reservoirs. The classical 
version of this problem was solved exactly by Bonetto et al ^33], using different methods. 
The quantum mechanical case was studied earlier by Vischer and Rich [32] and they 
obtained the thermal conductivity in the limit of weak coupling. The advantage of the 
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T 

Figure 2: Plot of the scaled thermal conductivity as a function of temperature (in units 
of hwc/ks) for u = 0.0. Inset shows the low temperature behaviour. 



present approach is that it's implementation is simple and straightforward. We obtain 
a general expression for the thermal conductivity, which in limiting cases gives both the 
classical result in [33] and the weak-coupling result in [32] • We also find that, at low 
temperatures, the temperature dependence of thermal conductivity in the case of finite 
coupling is completely different from the weak coupling case. In the classical case it was 
shown in [33] that the thermal conductivity of the harmonic chain with self-consistent 
reservoirs can also be obtained from the Kubo formula. An interesting problem would 
be to demonstrate this in the quantum-mechanical case. It is interesting to note that 
the self-consistent reservoirs are very similar to the Buttiker probes [33 [36] which have 
been used to model inelastic scatterering in electron transport. In the electron case they 
lead to Ohm's law being satisfied just as in the harmonic chain the introduction of self- 
consistent reservoirs leads to Fourier's law being satisfied. In fact we have recently shown 
how one can obtain Ohm's law using self-consistent reservoirs modeled microscopically by 
noninteracting electron baths ^37j. 

The Langevin equation approach has some advantages. For example, in the classical 
heat conduction case, it is easy to write Langevin equations for nonlinear systems and 
study them numerically. Also they might be useful in studying time dependent phenom- 
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ena. Examples of this are the treatment of quantum pumping in [3H] and the treatment 
of the question of approach to the non-equihbrium steady state in [22]. We feel that it is 
worthwhile to explore the possibility of using the quantum Langevin equation approach 
to the harder and more interesting problems involving interactions and time- dependent 
potentials in both the electron and phonon case. 
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A Green's function properties 

We will consider some properties of the Phonon Green's functions. We denote by G^{t) 
the full Green's function of the coupled system of wire and reservoirs. Let U and 
respectively denote the normal mode eigenvector and eigenvalue matrices satisfying the 
equations: 



U'^(!>U = fi^ U^MU = i . 



(A.l) 



We define the Green's function Q'^it) as 



sin [Qt) 



e{t) . 



(A.2) 



It satisfies the equation 



Mg+{t) + '^g+{t) = 5{t) I . 



(A.3) 



The Fourier transform Q^{uj) = dt ^+(t)e*'^* is thus given by 




(A.4) 



The isolated reservoir Green's functions are given by: 




1 




1 
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We can also represent G'^{uj) as follows: 



UtqUsQ 



[lu + UJQ + ie) {u -ujQ + ie) 



- E + ^ E UrQUsQ[5iu; - uq) + 5{u + ooq)] . (A.5) 

O Q O 



We will now express the wire-part of the full Green's function in terms of the uncoupled 
reservoir Green's functions. We write the equation for Q^{ijS) in the following form: 

-Mw {oj + iefi + Vl Vr 

Vl -Ml (u + ieyi + ^L 

-Mr (u + iefi + ^R 
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GwR \ 


( i 
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Grw 


Grl 


G^R J 


^ 








(A.6) 



From these equations we obtain the following expression for G^{u}): 

1 



-{u + ieY Mw + ^w-^t-^t ' 



where (a; 



= VlqUuj) Vl , 



(A.7) 



B Equilibrium properties 

In this section we will calculate the canonical ensemble expectation value oiK — X^X^ ) 
where the average is taken over the equilibrium density matrix of the entire coupled sys- 
tem of wire and reservoirs. Denoting by Zq the normal mode coordinates of the entire 
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system we get, for points i,j on the wire: 

= ^UiQUjQ{ZQ)eg 



Q 

TujJ 



Now from Eq. flA.7p we have: 

(G^)-^ - (G+ = (S+-S^) + (S+-S]j)+4^eu;MH. 

= 2 i {TL + TR)+4ieujMw 
^G+,-G^ = 2tG+,{TL + TR)G^ + Ateu;G:lyMwG^. (B.2) 

Hence we finally get: 

K!^ = r duj-[ (r^ + r^) G^ ^ coth( ^ ^ 



-co 



2A;bT^ 



+ / duj [G^MwG^ — coth(— -) . (B.3) 

J-oo 2 2/eBJ 

Since we eventually take the limit e — > 0, the second term is non- vanishing only when the 
equation 

Det[-uj'^Mw + ^w- ^ti^) - ^tH] = (B.4) 

has solutions for real u. These solutions correspond to the bound states [Ij of the coupled 
system of wire and reservoirs |22j. 
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